// Локальные
#include "linal/csrmatrix.h"

namespace linal
{
	//============================================================================================================================
	//============================================================================================================================
	void CSRMatrix::makeIdentity(int row)
	{
		m_di[row] = 1.0;

		for(int j = m_ig[row]; j < m_ig[row + 1]; j++)
			m_ggl[j] = 0;

                for(int j = m_ig2[row]; j < m_ig2[row + 1]; j++)
                        m_ggu[m_jg2[j]] = 0;

//                for(int j = 0; j < m_ig[m_Dimension]; j++)
//                        if(m_jg[j] == row)
//                        {
//                                m_ggu[j] = 0;
//                                int a = 5;
//                        }
	}

	//============================================================================================================================
	//============================================================================================================================
	CSRMatrix CSRMatrix::partFactorLU() const
	{
                CSRMatrix result(m_Dimension, m_ig, m_jg, m_ig2, m_jg2);
		for(int k = 0; k < m_Dimension; k++)
		{
			for(int j = m_ig[k]; j < m_ig[k + 1]; j++)
			{
				result.m_ggl[j] = m_ggl[j];
				result.m_ggu[j] = m_ggu[j];

				for(int i = m_ig[k]; i < j; i++)
					for(int m = m_ig[m_jg[j]]; m < m_ig[m_jg[j] + 1]; m++)
						if(m_jg[i] == m_jg[m])
						{
							result.m_ggl[j] -= result.m_ggl[i]*result.m_ggu[m];
							result.m_ggu[j] -= result.m_ggu[i]*result.m_ggl[m];
							m = m_ig[m_jg[j] + 1];
						}
						else
							if(m_jg[i] < m_jg[m])
								m = m_ig[m_jg[j] + 1];

				result.m_ggu[j] /= result.m_di[m_jg[j]];
			}

			result.m_di[k] = m_di[k];
			for(int i = m_ig[k]; i < m_ig[k + 1]; i++)
				result.m_di[k] -= result.m_ggl[i]*result.m_ggu[i];
		}
		return result;
	}

	//============================================================================================================================
	//============================================================================================================================
	double& CSRMatrix::operator()(int i, int j)
	{
		if(i == j)
			return m_di[i];
		else if(i > j)
		{
			int index = m_ig[i];
			while(m_jg[index] != j && index < m_ig[i + 1])
				index++;
			Q_ASSERT_X(index != m_ig[i + 1], "CSRMatrix::operator()", "The element hasn't been found");
			return m_ggl[index];
		}
		else
		{
			int index = m_ig[j];
			while(m_jg[index] != i && index < m_ig[j + 1])
				index++;
			Q_ASSERT_X(index != m_ig[j + 1], "CSRMatrix::operator()", "The element hasn't been found");
			return m_ggu[index];
		}
	}

	//============================================================================================================================
	//============================================================================================================================
	const double& CSRMatrix::operator()(int i, int j) const
	{
		if(i == j)
			return m_di[i];
		else if(i > j)
		{
			int index = m_ig[i];
			while(m_jg[index] != j && index < m_ig[i + 1])
				index++;
			Q_ASSERT_X(index != m_ig[i + 1], "CSRMatrix::operator()", "The element hasn't been found");
			return m_ggl[index];
		}
		else
		{
			int index = m_ig[j];
			while(m_jg[index] != i && index < m_ig[j + 1])
				index++;
			Q_ASSERT_X(index != m_ig[j + 1], "CSRMatrix::operator()", "The element hasn't been found");
			return m_ggu[index];
		}
	}

	//============================================================================================================================
	//============================================================================================================================
        CORE_EXPORT const Vector operator*(const CSRMatrix &matrix, const Vector &vector)
	{
		Q_ASSERT_X(matrix.dimension() == vector.dimension(), "operator*", "matrix and vector have different dimensions");

		Vector result(matrix.dimension());
		for(int i = 0; i < result.dimension(); i++)
			result[i] = vector[i]*matrix.di()[i];

		for(int i = 0; i < result.dimension(); i++)
			for(int j = matrix.ig()[i]; j < matrix.ig()[i + 1]; j++)
			{
				result[i] += vector[matrix.jg()[j]]*matrix.ggl()[j];
				result[matrix.jg()[j]] += vector[i]*matrix.ggu()[j];
			}

			return result;
	}
}
